Stability analysis of the D— dimensional nonlinear Schrodinger equation with trap and 

two- and three-body interactions 
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Considering the static solutions of the D— dimensional nonlinear Schrodinger equation with trap 
and attractive two-body interactions, the existence of stable solutions is limited to a maximum crit- 
ical number of particles, when D > 2. In case D — 2, we compare the variational approach with the 
exact numerical calculations. We show that, the addition of a positive three-body interaction allows 
stable solutions beyond the critical number. In this case, we also introduce a dynamical analysis of 
the conditions for the collapse. 
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Recent experiments on Bose Einstein Condensation 
(BEC) [U have brought great attention to its theoreti- 
cal formulation. Atomic traps are effectively described 
by the Ginzburg-Pitaevskii-Gross (GPG) formulation of 
the nonlinear Schrodinger equation (NLSE) M, which in- 
cludes two-body interaction. When the atoms have neg- 
ative two-body scattering lengths, a formula for the crit- 
ical maximum number of atoms was presented in ref . || . 
In ref. [^-^|, the formulation was extended in order to 
include the effective potential originated from the three- 
body interaction. In this case, in three-dimensions, it was 
shown that a kind of first order phase-transition occurs. 
In this connection, as also considered in the motivations 
given in , it is relevant to observe that recently it was 
reported the possibility of altering continuously the two- 
body scattering length, from positive to negative values, 
by means of an external magnetic field |7j] . Within such 
perspective, the two-body binding energy can be close 
to zero, and one can approach the so-called Efimov limit, 
which corresponds to an increasing number of three-body 
bound states ||. Near this limit, nontrivial consequences 
can occur in the dynamics of the condensate, such that 
one should also consider three-body effects in the effec- 
tive nonlinear potential. 

In the present work, we study the critical number of 
atoms in arbitrary D— dimensions, using a variational 
procedure; and also by an exact numerical approach in 
the case of dimension D = 2. The D— dimensional NLSE, 
with attractive two-body interactions, was previously an- 
alyzed in models of plasma and light waves in nonlinear 
media ||. The collapse conditions, in this case, were 
investigated without [|o| and with [jll| the harmonic po- 
tential term. In case of D =3, it was shown that a re- 
pulsive nonlinear three-body interaction term can extend 



considerably the critical limit for the existence of stable 
solutions 

Motivated by the observed high interest in stable so- 
lutions for arbitrary D, we look for variational solutions 
in a few significant cases (D =1,2,4 and 5) not previ- 
ously considered, when a three-body interaction term, 
parametrized by A3, is added to the effective non- linear 
interaction that contains a two-body attractive term. 
Our analysis also shows that, as in case of D = 3, a kind 
of first-order phase-transition can occur when D > 4, for 
certain cases of A3 > 0. In the present paper, we have 
also considered the approach given in |l2j, in order to 
study the stability conditions in the case of arbitrary D, 
when the non-linear interaction contains two (attractive) 
and three-body terms. 

In order to obtain an analytical approach and verify 
the validity of the variational Ritz method, we consider 
in detail the case of D = 2, with and without the three- 
body term, comparing the variational results with exact 
numerical calculations for some relevant physical observ- 
ables. In this case, we also discuss how the method given 
in p0[ | can be extended in order to approach analytically 
the exact value for the total energy. 

By extending the GPG formalism from three to D di- 
mensions, including two fl3]| and three-body interactions 
in the effective non-linear potential |^], we obtain 



in— — = 

dt 



2m 



A 2 M 2 + A 3 |V| 4 



(1) 



where ip = tp(f, t) is the wave- function normalized to the 
number of atoms N, to is the frequency of the trap har- 
monic potential and m is the mass of the atom. A2 and A3 
are, respectively, the strength of the two- and three-body 
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effective interaction, given in a D— dimensional space. such that 



\r\ is the hyperradius, such that r = 



D 



i and 

V = J2i=i &i~§F- i&i 1S the unit vector, with i = 1,2, ...D). 

The stationary solutions for the chemical potential fi 
are given by 



iU—— = aib. 
dt pv 



(2) 
5H 



Considering the general solution of eq.(lh, ifi— = 

^ dt dip* 

one can obtain the total energy E: 



E = / d D r H, with 



(3) 



W = |-|VV>| 2 
Zm 



H 2 + yH 4 + yM 6 . 



Here we consider only attractive two-body interaction, 
which is more interesting in the case of trapped atoms. 
For D =3, A2 = — 47rh 2 \a\/m, where a is the two-body 
scattering length and m is the mass of the atom. In the 
case of arbitrary D, A 2 has dimension of energy times 
L D , where L is a length scale in such space. However, a 
convenient redefinition of the wave-function in terms of 
dimensionless variables will absorb this constant, as will 
be shown. 

Our study will be concentrated on the ground state for 
a spherically symmetric potential. We first consider the 
case of A 3 = 0, using a variational procedure, with a trial 
Gaussian wave-function for ip(r), normalized to N, given 
by 



1pvar(r) = VN 



1 muj 



net 



D/4 



exp 



, (4) 



where a is a dimensionless variational parameter. From 
eq. (||), the corresponding expression for the total varia- 
tional energy can be expressed as 



w -ft ^ 



(5) 



D Da' 



4a 2 



G 3 v 3 n 2 



D 



4(2ir) D / 2 a D 6n D 3 D / 2 a 2D ' 

(6) 



where Qd is the solid angle in D dimensions, 
2ir D / 2 „ A 3 



G, 



r(D/2)' J "2(A 2 ) 2 
N 2A 2 fmuj\ D / 2 



frio 



flu huj ( fi ' 



(7) 



\(p(x)\ 2 d D x = N 



2|A 2 



(mui 



D/2 



(9) 



The dimensionless equation corresponding to eq. (0), can 
be rewritten as 



= 0, (10) 



where (3 = fj,/(hoj) is the dimensionless chemical poten- 
tial. From eqs. (g) and (||), the trial wave-function can 
be written as 



(f>var(x) = \/tMu ( — - ) exp ( 



D/4 



(11) 



The variational results, obtained by using the above 
expressions can be extended analytically to non-integer 
values of the dimension D. Minimization of the energy 
[eq. (||)], with respect to a 2 , is done numerically by 
sweeping over a 2 values. The results for the energy and 
the chemical potential are shown in Fig. 1. For each 
value of D, one can observe a critical number of atoms, 
N c , related to the critical parameter i/ c , only when D < 2. 
This critical limit corresponds to the cusps in the upper 
plot of Fig. 1 and is also observed using exact numerical 
calculation for D =3. It is also interesting to note that 
for D > 2 there are two branches of solutions for £ var and 
P, one stable and the other unstable. In the energy, the 
lower branch corresponds to stable solutions (minima), 
while the upper one gives unstable solutions (maxima). 




By using dimensionless variables, x = \J muj/h f, we re- 
define the wave-function ip as 



(S) 



FIG. 1. The variational energy £ (upper part) and chemi- 
cal potential j3 (lower part), as functions of the reduced num- 
ber of atoms v, for several values of the dimension D, indi- 
cated in each plot. All the quantities are in dimensionless 
units, as defined in the text. 
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The case with D = 2 is particularly interesting, as 
no unstable solutions exist and there are stable solutions 
only for v < 2, such that v c = 2. For D = 2, the mini- 
mization of eq. (||) in respect to a 2 leads to 



f — 



1 



(12) 



The behavior of v, and the corresponding critical limits, 
as one alters the dimension D, has other curious par- 
ticular results. For example, the critical limit u c has a 
minimum for D — 3 (vi^ > for all D). 

In conclusion of this part of our work, considering ar- 
bitrary D with A3 = 0, there are no stable solutions for 
eq. dJ), if the wave- function (f)(x), given by eq. (^|), is 
normalized to v > v c . Fig. 1 shows that this restriction 
is strongest for D — 3: v c is a minimum when compared 
with v c for D =/= 3. This is a relevant result, considering 
that v is directly proportional to the number of atoms. 
Also, it is observed that v c increases very fast for D > 3. 




FIG. 2. For D = 2, we present the variational and ex- 
act numerical calculations of the chemical potential total 
energy (£), mean-square-radius ((x 2 )), and central density 
(j<?!>(0)| 2 ), as a function of the reduced number of atoms v. All 
the quantities are in dimensionless units (see text). The solid 
line curves correspond to exact numerical results, while the 
dashed curves are the variational results. 

Next, we also solve equation ( |To| ) exactly employing 
the shooting and Runge-Kutta methods, and compare 
the results with the corresponding variational ones. In 
this case, we consider only the particular interesting case 
of D = 2, with A3 = 0. The results are shown in Fig. 2, 
for the chemical potential /?, the total energy £, mean- 
square-radius (a; 2 ), and the central density |^(0)| 2 . In or- 
der to numerically solve eq. (|l0|), in the s— wave, we first 



write it in terms of the single variable x = \J ix\ + x\) 
and consider the following boundary conditions: cf>'(0) = 
(where / stands for the derivative with respect to x) 
and <j>(x) — * Cexp(— a; 2 /2 + [f3 — 1] ln(x)) when x — > 00, 
where C is a constant to be determined. As observed in 
Fig. 2, the critical limit v c — 2 obtained analytically us- 
ing the variational approach should be compared with 
v c = 1.862, obtained by exact numerical calculation. 



This critical limit was first obtained by Weinstein |10|, in 
a non-linear approach with two-body term, without the 
trapping potential. The coincidence of the value with 
our exact calculation is due to the fact that at the criti- 
cal limit the mean square radius goes to zero. 
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FIG. 3. Exact numerical solutions for D =2, of the chem- 
ical potential, /3, in dimensionless units, given as function 
of the reduced number of atoms v, for different values of the 
three-body parameter G3, when the space dimension is D =2. 
As shown, only when G3 < the number of atoms is limited 
to certain critical number. 
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FIG. 4. The variational solutions for D =2, of the chemi- 
cal potential, corresponding to the exact results given in Fig. 
3. 
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We have also compared the results obtained by the 
variational approach with the exact numerical one, in the 
case of D — 2, for several values of the three-body inter- 
action term (positive and negative), as shown in Figs. 
3 and 4. In Fig. 3 we have the exact numerical ap- 
proach and in Fig. 4 we have the corresponding vari- 
ational results. By comparing the results we have for 
D — 2 (shown in Figs. 3 and 4) with the ones obtained 
in ref. || for D = 3, we should observe that no first order 
phase-transition exists in two dimensions. As observed in 
refs. for D = 3, a first-order phase-transition can oc- 
cur in trapped condensed states with negative two-body 
scattering length, when a repulsive three-body (quintic) 
term is added in the Hamiltonian. As shown in Figs. 
3 and 4, with G3 positive the range of stability for the 
number of atoms N can be increased indefinitely; with 
G3 negative this range is reduced. 

We can analyze the collapse conditions using "the virial 



theorem" approach |12|. The mean square radius, 



of a D— dimensional condensate, is given by 



d 2 (r 2 ) 
dt 2 



+ 4oj 2 (r 2 ) 



(13) 



1 

m 

where 



A(H)+X 2 (D-2)(\^ 



4A, 



(O) 



N 



d D r^{r,t)0^{f : t) 



(14) 



and (H) = E/N. When A3 = we obtain the equation 
derived in |12j| . 

We can also write the eq. ( |i"3| ) in dimensionless units, 
as it was done in eqs. (|^-|io|): 



dT 2 



Mx 2 ) = 



4£ 



where 
fir) 



A 2 D-2 
lAol 4 ' 



G. 



2/(r), 



D - 1 



(15) 



(16) 



Using the initial conditions for (x 2 ) and d(x 2 } /dr, where, 
for simplicity, we assume d(x 2 )/dT = 0, the solution of 
cq. (^||) is given by 



(x 2 ) = £ - 

V 



cos(2t) 
/(r')sin(2(r-r'))dr'. 



(17) 



The stability regions and the estimates for the collapse 
time can be obtained from the analysis of this solution, 
like as performed for the case A3 = in |H| . Let us ana- 
lyze the dynamics when D = 2. In this case, A2 does not 
appear explicitly in (x 2 ) and /(t) also does not depend 
on this parameter: 



1. For a positive G3, negative A2 and £ > we ob- 



serve that 



cannot be zero and the condensate 



is stable. The mean square radius of the conden- 
sate oscillates in time around a finite value. This is 
confirmed by the numerical simulations (see Figs. 
3 and 4). 

2. For a negative G3, positive A2 an analysis of stabil- 
ity like the one performed in ref. jlf| shows that 

a) When the total energy £ < 0, the condensate is 
unstable and the wavefields collapse in a finite time 
at any initial conditions; 

b) When £ > 0, as the function / (t) is negative, 
the contribution of the integral term for r < ir is 
negative. Then, we found the collapse condition as 



(x 2 )| >2 



£ 



(18) 



The same kind of analysis, for D > 2, is much involved 
in the present approach, as the sign of the function f(r) 
is not fixed at opposite signs for the parameters A2 and 

A3. " 

Some information about the dynamics of the collapse 
can also be obtaine d by u sing the techniques based on 

For instance, when D = 2, 



integral inequalities 1 10 



1( 



we can estimate the three-body term contribution in E, 
following the procedure given in |Kfl 



<r?#r < c 2 



2m 



d 2 r\^\ 2 ) = C 2 K 2 N, 
(19) 



where K is the kinetic energy and C2 is defined from the 
minimization of the functional 



J = 



(J d 2 r\V*jj\ 2 ) 2 (f d 2 r\t(>\ 2 ) 
fd 2 r\ifj\ e ' 



(20) 



Combining with the corresponding estimate for 
J d 2 f]ip\ 4 , we obtain E > E(K), where 

E(K) = K+ + —C\NK + ^-C 2 K 2 N. (21) 

v; ^ 4K 2 3 v; 

When A3 = we get the equation derived in p^] . Equa- 
tion ( |2l| ) should be compared with the corresponding 
variational expression ([5]), where the kinetic energy is 
given by K = Nhoj/(2a 2 ) and a is the width of the 
cloud. As we see, the expression for the energy (21) is 



very similar to the obtained by the variational approach. 
However, ( ^l| ) is valid for arbitrary time and describes 
the nonstationary dynamics. By using the variational 
expression (upper limit) for the ground-state, and the 
right-hand-side of eq. (|2l| ) (lower limit) , we can approach 
analytically the exact solution for the total energy 



E(K) < E < E v 



(22) 



4 



For a more deep insight to the problem of stability, we 
need to obtain the values of the constants C\ and G2. 
This problem requires a generalization of the method sug- 
gested by Weinstein in pOjl , to be considered in a future 
work. 
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FIG. 5. Variational solutions for the chemical potential (in 
dimensionless units) as functions of v, for D — 1 and different 
values of the three-body parameter G3. 



6 and 7, we present our results for the chemical potential 
as a function of v, for a set of given values of G3, in case 
of D = 1, 4 and 5. As one can observe in Fig. 5, even 
in case of D — 1 one can reach a critical maximum limit 
for v, when G3 is enough negative. For D — 4 and 5 
(Figs. 6 and 7), we observe similar picture of first-order 
phase-transition occurring for some specific values of G3. 
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FIG. 7. The same as in Fig. 5, for D = 5. 
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FIG. 6. The same as in Fig. 5, for D = 4. 

We should observe that exact numerical results, when 
G3 = 0, have already been considered in refs. H] (for 
D = 1 and D = 3), in g (for D = 3), and in @ (for 
D =2). In for D = 3, it was also considered the 

case with G3 7^ 0, and shown a kind of first-order phase- 
transition in the condensate. In the present work, we 
have extended the variational formalism, in case G3 =/= 0, 
for an arbitrary D— dimension. In the following Figs. 5, 



In conclusion, in the present work, we first studied 
the stability and the critical number of atoms in arbi- 
trary D— dimensions using a variational procedure, for 
the case we have two-body (attractive) and three-body 
contributions. This part extends a previous analysis done 
111 refs. @@. Next, we considered in more detail the 
case D — 2.. We compared the variational results with 
exact numerical calculations for the chemical potential, 
total energy, mean-square-radius and density. Finally, we 
extended numerically the approach for D = 2, including 
an effective three-body interaction term. We studied the 
sensibility of the critical numbers with respect to cor- 
rections in the non-linear interaction. The effective in- 
teraction considered in the equation contains a trapped 
harmonic interaction, and two nonlinear terms, propor- 
tional to the density \ip\ 2 (due to first-order two-body 
interaction) and to \ip\ 4 (due to first-order three-body in- 
teraction). We also verified, by a variational procedure, 
that a critical number of particles exists only for D > 2, 
when the nonlinear term of the NLSE contains just the 
cubic term. In case of D = 1 , a critical maximum number 
of atoms can exist with the addition of a negative quintic 
term (three-body) in the NLSE. In all cases where the 
number of atoms is limited, we observed that the addi- 
tion of a positive \ip\ 4 allows stable solutions beyond the 
critical number. We also introduced an analysis of the 
collapse conditions, using "the virial theorem" approach 
given in |lj] . The dynamics of the collapse was discussed 
in terms of the techniques developed in fj(| . In particu- 
lar, we showed how the exact energy can be approached 



5 



in the case of D 
tributions. 



2 with two and three-body term con- 
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